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ABSTRACT 

We present analytic radiative transfer solutions for the spectra of unresolved, spher- 
ically symmetric, centrally heated, dusty sources. We find that the dust thermal spec- 
trum possesses scaling relations that provide a natural classification for a broad range 
of sources, from low-mass protostars to dusty galaxies. In particular, we find that, 
given our assumptions, spectral energy distributions (SEDs) can be characterized by 
two distance- independent parameters, the luminosity-to-mass ratio, L/M, and the sur- 
face density, S, for a set of two functions, namely, the density profile and the opacity 
curve. The goal is to use SEDs as a diagnostic tool in inferring the large-scale phys- 
ical conditions in protostellar and extragalactic sources, and ultimately, evolutionary 
parameters. Our approach obviates the need to use SED templates in the millimeter 
to far-infrared region of the spectrum; this is a common practice in the extragalactic 
community that relies on observed correlations established at low redshift that may not 
necessarily extend to high redshift. Further, we demarcate the limited region of param- 
eter space in which density profiles can be inferred from the SED, which is of particular 
import in the protostellar community as competing theories of star formation are char- 
acterized by different density profiles. The functionality of our model is unique in that 
in provides for a self-consistent analytic solution that we have validated by comparison 
with a well-tested radiative transfer code (DUSTY) to find excellent agreement with nu- 
merical results over a parameter space that spans low-mass protostars to ultra-luminous 
infrared galaxies (ULIRGS). 

Subject headings: galaxies: formation — galaxies: starburst — infrared: galaxies — radiative 
transfer — stars: formation 



1. Introduction 

Stars and galaxies are born in dusty environments, shielded from view in the optical and often 
in the near-infrared. Radiation emitted by protostars and newly formed stars is absorbed by dust 



^ Department of Physics, University of California at Berkeley, Mail Code 7300, Berkeley, CA 94720 USA; 
sukanya® astro .berkeley.edu. 

^ Department of Astronomy, University of California at Berkeley, Mail Code 3411, Berkeley, CA 94720 USA; 
cmckee@astro.berkeley.edu. 



- 2 - 



and re-radiated at infrared and longer wavelengths. One of the primary tools for determining the 
physical parameters of these sources is the spectral energy distribution, or SED. 

The classic example of the use of SEDs to infer the nature of the underlying source is the 
classification of low-mass protostars based on the slope of their near-mid IR spectra (Lada & Wilking 
1984; Lada 1987; Adams, Lada, & Shu 1987): Class I objects, with d{vF^)/dhiv < in the 2-20 ^im. 
region of the spectrum, are identified with protostars; Class II objects, with < d{i'F,y)/d\'n.u < 2, 
are identified with classical T Tauri stars; and Class III objects, with 2 < d{vFy)/d\n.v < 3, are 
reddened stars approaching the main sequence. Subsequently, a fourth category. Class 0, was added 
for sources that are so embedded that it is generally not possible to measure the slope of the SED 
in the near- mid IR (Andre, Ward-Thompson, k, Barsony 1993). 

The physical interpretation of observations of massive protostars has proven to be a more 
challenging task. The task of interpretation is complicated by the greater extinction in high-mass 
star- forming regions and by the fact that massive stars often form in clusters. As a result, there 
is little consensus on cither the evolutionary parameters or the source parameters characteristic 
of massive star-forming regions. Even such a fundamental parameter as the formation time of a 
massive star has been uncertain by orders of magnitude. Nakano et al. (2000) interpreted near- 
IR spectroscopic data of a source in Orion as indicating a formation timescale of the order of 
10^ yr. Osorio, Lizano, & D'Alessio (1999, henceforth OLD99) estimated formation timescales of 
order 10^ yr from modeling the SEDs of hot cores. Estimates of formation timescales based on 
extrapolating the theory of low-mass star formation give formation times greater than 10^ yr, which 
is a significant fraction of the lifetime of the main-sequence lifetime of a massive star (McLaughlin 
& Pudritz 1997, Stahler, Palla, & Ho 2000). Subsequently, McKee & Tan (2002; 2003, henceforth, 
MT03) developed the Turbulent Core Model for high-mass star formation, which incorporates the 
effects of the supersonic turbulence and high pressures observed in massive star-forming regions 
(Plume et al. 1997) and predicts formation timescales of the order of 10^ yr. 

Other basic characteristics of high-mass star-forming regions, also remain uncertain. Only 
recently has it become possible to observe individual massive star-forming cores (Beuther &: Schilke 
2004, Cesaroni et al 1999, Fontani et al 2004), yet their properties remain uncertain because of the 
limited information on the SED. The density profile in high-mass star-forming regions is known 
only to be in the range of 1 to 2, with significant error bars (Shirley et al 2002, Beuther et al 2002, 
Jorgensen et al 2002, van der Tak et al 2000, Mueller et al 2002). On the other hand, OLD99 have 
argued for a more precise value, stating that their data on dust continuum emission from hot cores 
are best fit by so-called logatropic density profiles, with a power law index of 1. A self-consistent, 
analytic methodology for the inference of source and evolutionary parameters from the observed 
far-IR SEDs of protostellar regions would allow us to derive the large scale physical conditions of 
star-forming regions, and ultimately to discriminate between competing theories of star formation. 

A significant fraction of the star formation in the universe comes from dusty galaxies (Genzel 
k, Cesarsky 2000), including ultra-luminous infrared galaxies (ULIRGs), with luminosities between 
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8 jim and 1 mm in excess of lO^^L© (Soifer et al. 1984), and sub-millimeter galaxies at high 
redshifts (Smail et al 1997, Barger et al 1999a, Blain et al 2002), which appear to be even more 
luminous. SEDs are a primary tool for inferring the properties of these sources. The lack of data 
on high-redshift sources has led to the use of SED templates from low-redshift sources to infer star 
formation rates in high-redshift sources (e.g., Xu et al 2001, Carilli &; Yun 1999, Yun Sz Carilli 
2002, henceforth YC02), but it is not known whether these templates are valid. The discovery 
of massive, luminous sub-millimeter galaxies may warrant a revision of the star formation history 
of galaxies at high redshift (Blain et al. 2002) and has been suggested as being problematic for 
hierarchical formation scenarios (Genzel et al. 2003). The importance of understanding the SEDs 
of high-redshift galaxies is highlighted by the work of Chapman et al. (2004), who have proposed 
that higher temperatures in these sources mean that current sub-mm surveys may have missed 
more than half of the most luminous, dusty galaxies at z ~ 2. 

The simplest model for an SED is to assume that it is due to an isothermal distribution of 
optically thin dust (Hildebrand 1983). The next level of sophistication is to allow for the dust to 
be optically thick above some critical frequency, while still considering the entire dust envelope 
as being characterized by a single temperature, e.g., Yun & Carilli 2002, henceforth YC02. A 
difficulty with such single-temperature models is that they often require /? to be less than the value 
appropriate for uncoagulated grain models (/? ~ 2; Weingartner & Draine 2001). This difficulty 
can be overcome if the spectrum is assumed to result from the superposition of two blackbodies at 
different temperatures (Dunne &; Eales 2001). 

In reality, the dust temperature varies continuously. Variation of the dust temperature has been 
taken into account in work on low- mass star formation beginning with the pioneering work of Larson 
(1969), who introduced a physically motivated approximation for the temperature profile. Adams & 
Shu (1985) presented an approximate numerical radiative transfer model based on this form of the 
temperature profile and showed that they could approximately satisfy radiative equilibrium. With 
this model, they inferred stellar masses and accretion rates for their favored collapse model, the 
singular isothermal sphere. This paper heralded the beginning of a stream of papers on more refined 
numerical modeling of low-mass protostellar spectra, with inclusion of effects for more evolved 
sources. Kenyon, Calvet, & Hartmann (1993) developed an approximate numerical approach in 
which the temperature profile is calculated from the radiative diffusion equation in the optically 
thick part of the envelope and from radiative equilibrium in the optically thin part. OLD99 adopted 
this approach in their calculation of the SEDs of massive protostars. Numerical modeling has also 
enabled the modeling of axisymmetric sources, including disks (Efstathiou & Rowan-Robinson 
1991; Kenyon et al 1993; Whitney et al 2003). However, analytic radiative transfer models for even 
the simplest spherically symmetric systems remain rare. An exception is that of Adams (1991), 
who presented an analytic solution for the specific intensity of protostellar cores at millimeter and 
sub-millimeter wavelengths; however, he did not evaluate the accuracy of his method. 

While numerical models permit one to calculate the SED of a given source with exquisite 
accuracy, what is lacking is any general understanding of how the SED depends on the underlying 
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source parameters. A significant step in addressing this problem was taken by Ivezic &; Elitzur (1997, 
henceforth IE97), who emphasized the importance of scahng relations in determining the spectrum 
of dusty sources. They showed that the spectrum of a spherically symmetric, dusty source is 
determined by four parameters: the dust destruction radius, Rdd, which depends on the luminosity 
of the source (see Appendix A); the thickness of the shell around the source, y = Rc/Rdd^ where 
Rc is the radius of the cloud in which the source is embedded; the optical depth through the shell 
at some frequency: and the temperature of the central source. In addition, the spectrum depends 
on two functions, the density distribution, p(r), and the opacity, k^- Using these scaling properties 
of dust emission, they developed a numerical code, DUSTY, that is very useful in inferring the 
physical conditions in dusty sources. 

Our objective in this paper is to develop an analytic theory for the SED of spherically symmet- 
ric, dusty sources. We assume that the dusty envelope surrounding the central source of radiation 
is sufficiently opaque that the resultant SED is approximately independent of the temperature of 
the central source and of the properties of the dust destruction front. Within our parameter space, 
lOOg cm~^ ^ ^ ^ 0.03g cm~^, the optical depth is large enough that essentially all the stellar radi- 
ation is absorbed by dust and re-emittcd. Because of the large optical depth, the near-infrared and 
optical spectrum, which does depend on the source spectrum, is heavily attenuated. For a given 
form for the density distribution (e.g., a power-law) and the opacity (e.g., that of Weingartner & 
Draine 2001), the emergent spectrum depends on three parameters: the luminosity, L; the dust 
mass, M(just = -^/■^dust) where .^dust is the mass fraction of dust; and the radius of the source, 
Rc- The shape of the SED is independent of the distance, and therefore depends on only two 
parameters, which we take to be the luminosity-to-masss ratio, L/M, and the surface density of 
the source, S = M/{ttR'^). As we shall show, it is possible to infer these two distance-independent 
parameters, and therefore the complete shape of the SED, from just two colors (provided, of course, 
that the assumptions underlying our model are correct). Since the pressure in a self-gravitating 
gas is of order GS^, determination of the surface density allows one to infer the pressure in the 
source. If the distance is also known, we can infer the luminosity, dust mass, and physical size of a 
protostcllar region-even if the source is unresolved. However, it is generally not feasible to infer the 
density profile from the far-IR SED of an unresolved source; as we shall see later, this is feasible 
for extended envelopes, which are large compared to the radius of the effective photosphere, and 
for envelopes that emit most of their radiation at wavelengths shorter than 30/Ltm. 

We find that the spectra are characterized by three frequency regimes: Low frequencies, which 
are always in the Rayleigh- Jeans limit [hu < kT{Rc)]; intermediate frequencies, which are not 
necessarily in the Rayleigh-Jeans limit, but are low enough that the envelope is transparent; and 
high frequencies, where the spectrum is determined by competition between opacity effects and 
the Wien cutoff. Having found formal solutions in these three regimes, we present a joint solution 
within a conceptually simple and physically motivated framework in which the emission in each 
frequency regime comes from a shell of some thickness, centered at some radius, and is attenuated 
by the intervening optical depth. We adopt a self-consistent temperature profile that characterizes 
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the emission coming from the vicinity of an effective Rosseland photosphere. The slope of this 
profile is set by the condition that the emergent luminosity equal the input luminosity. We have 
tested the accuracy of our analytic solution by making detailed comparisons with DUSTY. 

The organization of the paper is as follows: in §2 wc outline the general formulation of the 
problem, and specify the forms of the density profiles that we consider as well as the dust opacity. 
In §3 we introduce the characteristic parameters that are an integral part of our formalism. In §4 
we present the formal solution to the problem in the three frequency regimes that are characteristic 
of the dust thermal spectrum, and in §5 we present the joint solution for the SED. In §6 we 
discuss luminosity conservation and the resultant form of the temperature profile. §7 presents 
the accuracy and range of applicability of our solution. §8 is dedicated to an explanation of 
our results and analysis, with particular attention to the shape of the SED and its dependence 
on two parameters, the feasibility of inferring density profiles from the SED, and the emergent 
three-component spectrum for highly extended envelopes. In Paper II, we shall present the far-IR 
SEDs and inferred source and evolutionary parameters for a broad range of sources, from low mass 
protostars, to massive protostars and ULIRGS. 



2. Formulation of the Problem 

We consider a centrally concentrated source of radiation surrounded by a homogenous, spher- 
ical distribution of dust. Consideration of a central source of radiation, while appropriate for 
protostcllar sources, is approximately valid for ULIRGs, particularly if the far-IR emission is pre- 
dominantly powered by an extended starburst. As such, our method is more applicable for ULIRGs 
largely powered by dust-enshrouded AGN or compact star clusters, and for super-star clusters in 
starburst galaxies and ULIRGs. A relevant finding here is that of Soifer (et al 2000) - they find that 
that a large fraction of the mid-infrared emission stems from very compact systems. If the source 
is unresolved, our method is applicable when the dust temperature in the beam is dominated by 
the central source. If there are multiple sources present, or if the background temperature is signif- 
icant, then the angular resolution needs to be sufficient to meet the above condition. If the angular 
resolution of the beam is sufficient to resolve the source, our method can be applied to compute the 
emergent SED given that all the flux from the source is included in the beam. Our assumption of 
spherical symmetry limits us to consideration of young protostellar sources in which any accretion 
disk is small enough with respect to the surrounding envelope that it does not significantly affect 
the SED. We approximate the density distribution by a power law in radius. Wc also assume that 
the dust is homogeneous; the effect of clumping will be considered in a future paper (Chakrabarti 
& McKee, Paper III). 

We neglect scattering of radiation by dust grains and consider only the thermal emission in 
computing the emergent SED, as the scattering efficiency scales as X"^ and is thus important only 
at very short wavelengths. We consider dust shells sufficiently opaque that the dust destruction 
front is highly obscured, and thus does not significantly affect the far-IR spectrum (see Appendix 
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A) . This also implies that the observed SED is independent of the spectrum of the central source 

of radiation. We consider the opacity to be a function of frequency only; we do not consider the 
temperature dependence of the opacity that can be caused by evaporation of grains. Although we 
do not treat the temperature dependence of the opacity explicitly, we give a simple prescription to 
implement it in Paper II. These assumptions are discussed in more detail in §7. 

We find that we get excellent agreement with the numerical solution from DUSTY for the 
far-IR SED with the adoption of a power law temperature profile. We have specified the variation 
of the slope of the temperature profile over a parameter space that spans a range of optical depths 
of a factor of ~ 1000, through a combination of heuristic arguments and numerical calibration with 
DUSTY. Thus, although the actual temperature profile is not a pure power law, the mm to far-IR 
emission can be well described by an effective power-law temperature profile that characterizes 
the emission coming from close to the Rosseland photosphere. Most of the observed emission 
originates outside the Rosseland photosphere, since emission at high frequencies is attenuated by 
the intervening optical depth, and emission at low frequencies comes from outside the r = 1 surface 
at the peak of the SED. 

As explained in §6, we adopt a temperature profile that is a result of imposing the self- 
consistency criterion that the input luminosity be equal to the emergent luminosity. We find that, 
for a given density profile and dust model, the slope of the temperature profile is a function of one 
dimensionless parameter, insofar as most of the emitted flux is longwards of A ~ 30 /xm, where the 
opacity is approximately a power law in frequency. 

2.1. Density Profile 

Observations of low and high mass star forming regions (Shirley et al 2000, Shirley et al 
2002, Beuther et al 2002, Jorgensen et al 2000, van der Tak et al 2000, Mueller et al 2002) have 
found power law density profiles, p{r) oc r'^f, with density power law index, kp, in the range 
of 1 to 2. As discussed in more detail in §8, within our formalism we approximate kp = 1 with 
kp = 1.1, as the emitted spectra agree nearly exactly. We note that most of these authors have 
found power law density profiles from a combination of fitting to the observed SED and intensity 
profiles. (We defer detailed discussion of the feasibility of inferring density profiles from the SED 
alone to §8.) The observed prevalence of power law density profiles in protostellar envelopes can 
be heuristically understood in the context of a self-similar, turbulent medium (MT03). The self- 
similarity should be broken only on small scales, by thermal motions (MT03), as they are observed 
to do in regions of low-mass star formation (Blitz & Williams 1997). However, such a picture of 
a time-stationary, spherical, self-gravitating turbulent structure is necessarily highly approximate, 
and large fluctuations, both spatial and temporal, are to be expected. We defer consideration of 
these fluctuations, which correspond to dumpiness, to a future paper. 

It is important to note that the observations of massive star forming regions cited above probed 
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scales of clumps, which are massive enough to form a cluster of stars, rather than of cores, which 
will form a single star or binary. (We follow the clump-core terminology introduced in Williams, 
Blitz &: McKee 2000.) However, as we shall show in Paper II, the photosphere of a high-mass 
protostar is within its natal core. Cores are too small to have measured density profiles yet. We 
assume that they also have power-law density profiles, which is consistent with MTOS's assumption 
that clumps and cores are part of a self-similar structure. We assume that the cores have a definite 
outer radius at Rc, beyond which the density drops rapidly; we assume that the emission from 
beyond Rc is negligible. 

While the simplest SED models of star-forming galaxies do not consider a density variation 
(e.g. YC02), more sophisticated radiative transfer models (Efstathiou et al 2000) have considered 
starburst galaxies as an ensemble of optically thick clouds heated by the newly formed stars. The 
assumption of power-law density profiles for ULIRGs is a very approximate representation of the 
complex morphology of merging systems; it may be a good approximation for the clouds that 
comprise the ULIRGs, however. 

2.2. Dust Opacity 

The results presented in this paper are based on the Weingartner & Draine (2001) (henceforth 
WDOl) Rv = 5.5 dust model. WDOl is an extension of the original Draine &; Lee (1984) dust 

model, with a size distribution developed to reproduce the observed extinction curve for a variety 
of environments, parametrized by the ratio of the visual extinction to reddening, Ry. For the 
diffuse ISM, Ry is observed to be approximately 3. Higher values have been observed for dense 
clouds (Strafella et al 2001, Kandori et al 2003, Vrba et al 1993) and star-forming galaxies (Calzetti 
2000). 

The WDOl model has a simple composition, consisting of carbonaceous grains and silicates. 
WDOl's Rv = 5.5 model has a substantial depletion of the smallest carbanaceous grains relative to 
the Rv = 3.1 model, leading to a difference in the opacity curves in the near-IR. However, WDOl's 
best fit models for Ry = 5.5 and -Ry=3.1 have no diff'crence in the far-IR extinction curves, and 
thus imply no substantial grain growth for sizes on the large end of the size distribution, i.e., 
a ~ 0.1/Ltm. As we later show, for densities typical of protostellar regions and dusty galaxies, it is 
unlikely that significant coagulation occurs within a few free-fall times to affect the spectrum in the 
mm to far-IR region of the spectrum. For uncoagulated grain models such as WDOl, the mm to 
far-IR variation of the opacity with frequency (in the range of 3 mm to 30 /xm) is well represented 
as a power law with slope, f3 = 2. 

The opacity normalization per gram of dust, k^q, depends on the metallicity. Let be the 
mass fraction of dust, equal to 1/105.1 for solar abundances (WDOl) and 5 be the dust-to-gas 
ratio relative to solar. At a fiducial wavelength of Aq = lOO/xm, WDOl's opacity is k^/q = 0.275, 
independent of Ry. The WDOl models reproduce the observed extinction curves for the Milky 
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Way, as well as regions of low metallicity, such as the LMC and the SMC. Since WDOl is the 
simplest grain model that is able to reproduce the observed observed extinction over a wide range 
of metallicities, we have adopted it as our fiducial dust model. 

The principal omission in the WDOl model is the lack of ice mantles. Spectroscopic studies 
have identified ice spectral features in protostcllar environments and in dusty galaxies (Allamandola 
et al 1992, Tielens et al 1984, Spoon et al 2002). Based on observations with the Submillimeter 
Wave Astronomy Satellite, Bergin et al. (2000) have argued that most of the oxygen in star-forming 
regions is frozen out onto the dust grains in the form of molecular ices. Preibisch et al. (1993) 
have modeled dirty ice mantles, and find significant variations in the opacity, depending on how the 
mantles are distributed among the grains and on the volume fraction of carbonaceous material; a 
typical increase in the far-infrared opacity is a factor ^ 2. Pollack et al (1994) considered detailed 
grain compositions, including the contribution from water ice. If most of the oxygen is in the form 
of water ice, they find that Ki,q 2± 1 cm^ g"-*^ at 100 /Ltm; however, at = Aq = 1 mm, the opacity 
for spherical grains of average radius 1 /xm is 4.2 x 10~^ cm^ g~^, not that different from WDOl 
(k^,,) = 3.0 X 10^'^ cm62 g""*^). However, the steep power law index for the opacity, (3 = 2.7, found 
by Pollack et al. docs not appear to be consistent with observations. Pollack et al. show that 
for typical conditions in regions of high-mass star formation (nn ^ 10^ cm"^), ice sublimates at 
T 2± 110 K. We conclude that sigificant uncertainties remain in the effect of ice mantles on the 
far-infrared opacity. We show later in §3 that the effective increase in opacity normalization that 
ice mantles produce may be simply treated within our formalism. 

A further complication to consider is the possible growth of grains due to coagulation (Pollack 
et al. 1994; Ossenkopf &; Henning 1994). The time scale for coagulation depends on the mean 
relative speed of the grains, which is likely due by turbulence (Draine 1985; Lazarian & Yan 2002) 

in dense cores. The outcome of a collision between grains depends on whether or not the grains 
stick together, i.e., if the surface potential energy is comparable to or larger than the kinetic energy. 
Recent numerical experiments (Poppe et al 2000) have shown that relative velocities of order Im/s 
allow sticking probability of order unity, with a sharp decline at v ~ 5m/s. Therefore, the timescale 
for grain coagulation is approximately 



where we have considered the approximate case of coagulation of identical, spherical grains of size 
a. To determine whether coagulation will be effective within timescales relevant for protostcllar 
evolution, we compare this coagulation time scale to the free fall time. For coagulation to affect 
the spectrum in the far-IR range would require grain growth to the Rayleigh limit, i.e., to sizes 
a ~ A/27r, a > 5 fim, (to affect emission at A > 30jUm). This gives: 



Therefore, for typical estimates of the densities, nn ~ 10^ cm ^, we see that the coagulation 
timescale is more than several hundred times the free-fall time for a ~ 5/xm, which corresponds to 




(1) 





-9- 



A = 30/Ltm/27r. This conclusion is consistent with the results of Chokshi et al. (1993), who found 

a relatively small increase in grain size (approximately a factor of 2) in dense cores. We conclude 
that coagulation in star-forming regions is unlikely to significantly alter the far-infrared opacity. 



3. Characteristic Parameters 



Stars have well-defined photospheres, but dust clouds do not. Nonetheless, the SED of a dusty 
source can be described in terms of characteristic radius, Rch, and characteristic temperature, T^h, 
such that 

L = ATTLRlaT^^ , (3) 

where L is a number of order unity determined below in order to secure better agreement with the 
actual SED. We determine i?ch and T^h by requiring that a characteristic optical depth at frequency 
fch = kTch/h is unity, 

Teh = Kch J p{r)dr = 7 = 1 ; (4) 



kp-l 

note this characteristic optical depth ignores the cut off in the density at the edge of the cloud, 
which is at Rc- Rch and Teh are therefore the approximate photospheric radius and temperature; 
more accurate photospheric values are given in §7 below. 

We now express the characteristic parameters in terms of the physical source parameters, the 
surface mass density, S = M/nR^, and the luminosity to mass ratio, L/M. We assume that u^h 
lies within the frequency regime in which the opacity is a power law. 



1^1^ = Ki^oC^/^o)^ (30 /Ltm < A < 1 mm). 



(5) 



Solving equations (3) and (4), we evaluate the two parameters that we use to describe our solutions, 



Rr. = 



Rc 
Rch 



(L/M)S(4+/3)//5 



4aL 



(3 



4(/c, - l)T^ 



i/3 



4//3 ^ 2/3+4(fep-l) 



(6) 



and 



ch 



j L/M 
I 



4(fcp - l)n 
(3 kp)ni,Q 



2/3+4(fep-l) 



(7) 



In general, L is a function of Rc, the form of which is specified in §6. It is important to note that 
while we have defined the characteristic parameters in terms of a power law opacity for simplicity, 
our solution for the emergent SED is valid for an arbitrary opacity curve, as we show later. 

The utility of equations (6) and (7) is to allow an analytic translation between the SED 
variables, Rc and T^h, which govern the shape of the SED, and the source parameters, L/M and S. 
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The luminosity-to-mass ratio and the surface density are independent of distance, and therefore, 
the SED variables, Rc and T^h, are also. Once the SED variables are determined from two observed 
color ratios, we may then immediately solve for L/M and S using these equations. 

The SED variable Rc = Rc/Rch is a dimensionless measure of the size of the region that the 
observed far-IR SED probes. When it is large, the observed emission arises from a wide range of 
radii; when it is somewhat larger than unity, the emission comes from a narrow range. Values of 
Rc less than unity are generally not meaningful since the resulting spectrum is sensitive to our 
assumption that the emission aX r > Rc drops to zero. As we shall see below, dense dust shells, 
such as those characteristic of massive protostars, starbursts and AGN, have smaller values of 
Rc than lower density dust shells, such as those characteristic of low-mass protostellar envelopes. 
Thus, for low-mass protostellar envelopes, we probe the source function over a larger region than 
for high-mass envelopes. 

The angular size of the photosphere is about 

^ch = , (8) 

where D is the distance to the source. If the total flux, F = L/AttD^, is known, then the angular 
size is 



from equation (3). Thus, since Teh determines the surface brightness of the source, it is possible 
to infer the angular size of the source even if it is unresolved. If the distance to the source is also 
known, we may then further solve for the size of the source, as well as the mass and luminosity. 

In Figure la, we have depicted the typical locations of various sources on the L/M vs S plot. 
Isolated low-mass protostars appear typically on the lower left corner of this diagram, i.e., at low 
surface density, S ~ 0.05 g cm~^, and low luminosity to mass ratios, L/M ~ 1 Lq/Mq (Jorgensen 
et al 2002). In clusters, E is higher since the pressures are higher (MT03). High-mass protostars, 
ULIRGs and super-star clusters are typically in regions of higher surface density, E ^ 1 g/cm^ 
(Plume et al 1997, MT03), and thus appear on the right-hand side of this diagram. ULIRGs and 
super-star clusters have comparable L/M 10 Lq/Mq on average (Downes & Solomon 1998), 
(Gilbert &: Graham 2002, Gilbert 2002). Two sets of lines of constant T^h and Rc are overlaid 
on this diagram for our fiducial density profile, kp = 3/2, and our adopted dust model WDOl. 
As discussed above, ice-coated grains produce an effective change in the opacity normalization 
of a factor of ~ 2. Our scaling relations (6) and (7) show that this leads to a small change in 
the governing SED variables. The difference relative to WDOl's normalization when ice-coated 
grains are considered is depicted in Figure lb for our fiducial density profile. A more quantitative 
discussion of the characteristic parameters is relegated to §7. Figures 2a, b, and c present fits to 
observed data for a low, high-mass protostar and ULIRG, respectively. Paper II gives a detailed 
explanation of the application of our methodology to star-forming system, along with SED fits to 
about a dozen sources, from low-mass protostars to ULIRGs. 




(9) 
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4. Evaluation of Li, 

We compute the spectral luminosity as 

= 47r j,Mr)4Trr^dr , (10) 



where jjy is the emissivity and fuir) is the escape probability, which is evaluated in Appendix B. The 
inner boundary of the emitting shell is at Rddj the dust destruction radius, and the outer boundary 
is at Rc, the radius of the cloud. At high densities, the dust and gas temperatures are the same, 
but at low densities they may differ; in our equations, T always refers to the dust temperature. 
We assume that the dust emissivity is given by the LTE expression, ji, = K,i,p{r)Bi,(T). This 
approximation breaks down at the dust destruction front, where the absorption of UV photons can 
lead to transient heating of small grains (Drainc &: Anderson 1985). However, since we assume 
that the cloud is sufficiently opaque that the dust destruction front is shielded from view, the LTE 
approximation should be valid. 

In order to evaluate the luminosity, we must specify the temperature profile. We assume that 
it can be approximated as a power law in the vicinity of the photosphere, 

T = reh(^-^) =T,^f-^^ , (11) 

where f = r/Rch. We determine kx to zeroth-ordcr by imposing the self-consistency condition that 
the input luminosity equal the emergent luminosity - the procedure is described in detail in §6. Note 
that this assumed form for the temperature carries the implicit assumption that the temperature 
at Rch is Teh', we choose L{Rc) choose so as to improve the accuracy of this approximation. Written 
in dimensionless notation, our expression for the luminosity then becomes 

= 47ri?,\47r (^^) ^.^'(fcp - 1)/ , (12) 



where 



1= / ■'7.:, ^ -dr , (13) 

exp(z/r*^T) _ 1 ^ ^ 

9 = i^/t'ch) Ki/ = Kjy/Kjy^j^, and hu/kT^r) = vf^'^ . Note that since the integral is performed over 
position, we may take the opacity term, k^,, outside the integral since we have assumed the opacity 
is not a function of position. 

In order to obtain an approximate analytic evaluation of /, we have found it necessary to 
consider three distinct frequency regimes, which we denote as low, intermediate, and high. The 
low and intermediate frequencies are optically thin. Low frequencies are in the Rayleigh-Jeans 
portion of the spectrum throughout the envelope [hv < kT(Rc)]. The low-frequency emission comes 
predominantly from the outer parts of the shell, as it is proportional to the mass. Intermediate 
frequencies are in the Wien portion of the spectrum in the outer envelope, although not near the 
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photosphere. The temperature dependence of the intermediate frequency region causes the higher 
components of this frequency regime to emanate from deeper in the envelope, in a sense specified 
by the temperature variation of the envelope. High-frcqucncics are in the Wien portion of the 
spectrum at the photosphere; the emission originates from a location that is due to a tug-of-war 
between the hotter temperature in the interior and the intervening optical depth that has to be 
traversed. 

Once we obtain an approximate expression for in each of these frequency regimes, it is 
necessary to knit them together into a single expression. To do this, we have found it convenient to 
introduce a "shell" model, in which the emission at each frequency comes from a shell of thickness 
Arm(z/) centered at a radius r^(z^), with a source function exp [—hvlkT{fjn)] located at 

an optical depth T,^{rrn)- 

= 47rRli,4n (^^^) K^i>Hkp - l)?^^"^" Af^ exp 

where the optical depth from r to the surface of the cloud is 

T. = (f-'"-^' - Rc'"^') . (15) 

We now proceed to evaluate the parameters of this intuitively transparent form for the spectral 
luminosity. 



hu 



(14) 



4.1. Low and Intermediate Frequency Regimes (rjy <C 1) 

The low and intermediate frequency regimes are characterized by Tj^ <C 1, so that the escape 
fraction is approximately unity, /jy(r) ~ 1. Low frequencies are in the Ray leigh- Jeans part of the 
spectrum throughout the cloud, so that hv/kT = i'f^'^ <C 1 everywhere and the exponential in / 
can be expanded, 

= {3-\-kTy ^''^ 

where we have assumed that i?dd^~'^''~^'^^ ^ Rc so that the dust destruction radius does not 

significantly affect the spectrum; equivalently, the emission is dominated by the outer part of the 
envelope. We note that this requires kr < S — kp. 

Intermediate frequencies are in the Wien part of the spectrum in the outer part of the cloud 
[hu /kT{Rc) ^ 1), and as a result the upper limit integration can be set to infinity. Since we have 
assumed that the dust destruction front does not affect the spectrum, we then find 



ZJIA c f \ (17) 



(Gradshteyn & Ryzhik p. 349). 
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Now we combine the formal expressions for the low and intermediate frequency fluxes into a 
single expression by observing the scaling behavior of each. In this case, we require a harmonic 
mean to recover each limiting case, = -^iow"'"-^int • "^^^ resulting expression for the luminosity 

is 



i^,.,iow-int - leir'^ikp - 1) ( ) «<^^chZ>^ 



rc 



i3-kp- kT)i>TCRc 



-kr) 



(18) 

where the argument {?> — kp)/kT for the Gamma and Zeta functions has been suppressed for clarity. 
This prescription holds for kp in the range of 1 to 2 that we are considering. 

Recall that in our shell model (eq. 14), we characterize the emission at each frequency as 
coming from an optically thin shell centered at a radius r^(z^) and having width ArTO(i')- To apply 
this equation, wc must dotGriniiiG tli6 shell radius cind. tliickiiBSS. Tli6 sIibII ra,dius 

is also known 

as the contribution function (e.g., Gray p. 278). As discussed below equation (16), most of the 
low- frequency emission comes from the outer part of the shell, so we set r^, low = Rc- To determine 
the shell radius at intermediate frequencies, we assume that hv /kT{fm,\nt) is a constant, 



hu 



kTif. 



m, int ) 



m, int 



c, 



(19) 



where C is determined from numerical calibration in §7 and §8.1. To recover the limiting cases, 
the shell radius for both low and intermediate frequencies can be approximated by the harmonic 
mean, so that 

rm, low-int - ^^~l/kr+C^/kT ' 
As we discuss in detail in §8.1, when Tm^iow — '"m, int; 

the slope of the emitted spectrum changes 
from the standard Rayleigh-Jeans slope to the flatter intermediate frequency slope. This break 
frequency, i.e., the frequency at which the character of the emission changes from being dominated 
by the cool material on the outside to the emission coming from deeper into the envelope in a sense 
specified by the temperature gradient, is then given by: 



^reak 



CR. 



In terms of the temperature at the outer edge of the envelope, the break frequency is: 

kT{Rc) 



^reak 



h 



(21) 



(22) 



In order to recover equation (18) for the luminosity at low and intermediate frequencies, we 
must set the thickness of this emitting region to be 



m, low — int) 



(S-kp- kT)i>rCRi 



-{3—kp—kT) 



1 



p2 kp 
m, low— int 



(23) 
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4.2. High-Frequency Regime 

In the high-frequency regime, the optical depth can be larger than unity and it is necessary 
to develop an escape fraction appropriate for spherical geometry. In Appendix B we show that an 
appropriate escape fraction for spherical geometry is: 

" l + 2.„(|^)(f-'=--PiJ-.-) • 

where t^, is the radial optical depth and is given in equation (15). We note that for a given dust 
model and density profile, the optical depth depends on Tch in a scale-free manner so long as we 
consider wavelengths longwards of 30 //m, where Hp oc z/^ is a power-law in frequency. 

We compute the high-frequency luminosity by the method of steepest descent. This method 
is particularly useful when the integrand is composed of a rapidly varying function and a slowly 
varying function, such as an exponential term multiplied by an algebraic term, which is indeed 
the structure of the high-frequency integrand. The method of steepest descent approximates the 
integral as the algebraic term evaluated at the local maximum of the exponential and a gaussian 
centered at the local maximum, having a width that is proportional to the second derivative of 
the argument of the exponential. The narrowness of this gaussian is an indication of how local the 
contribution to the integrand really is. The narrower the gaussian, the better the approximation. 

Inserting our expression for the escape fraction into equations (12) and (13) we find 

where the argument of the exponential term is 

h{r) = vf^^ + (r^-'''' - kl''""^ . (26) 

The location of the maximum of the exponential term [i.e., the minimum of h(r)\ gives the 
shell radius at high frequencies according to the method of steepest descent. 



'"m, high, steep 



f^vikp 1) 



l/(kT+kp-l) 

(27) 



Physically, this means that the high-frequency photons come from a location in the shell that is due 
to a competition between the optical depth and the temperature gradient, with the temperature 
gradient driving fra{v) inwards, while the optical depth drives fm{i^) outwards. 

In principle, this expression for the shell radius can give r^high > whereas in fact fm Rc 
as the shell envelope becomes very opaque. Thus, in general the shell radius is given by 

^m, high = Min(i?c, Vm, high, steep) • (28) 
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It is important to note that, in finding the contribution function for the high-frequency regime, 

we made no assumption about the frequency dependence of the opacity curve, and thus wc may 
use a tabulated, reaUstic opacity curve to find for the high-frequency regime. We have assumed 
that Ky is independent of position for simplicity, but in principle equation (27) could be generalized 
to include temperature dependent opacity, At,^ [T(r)]. 

Evaluating the integral in equation (25) by the method of steepest descent, we find 

+ ^'^i^ (^fcp+l j V m,high ''m,high-"c ) 

where hm = h{fm, high) and 

C = i>kT{kT - l)rt-j:Jh + KiK - l)'^. V^fgi • (30) 

We can express the luminosity in the form of the shell model (eq. 14) if we identify the thickness 
as 

i'^^W' . 

^?-m,high - ^ ^ ^, (fc^_i) _ ,2 p-'tp-lV ^ ' 

1 + ^Ki/ j^^j^-^ high ^m, high-^-c j 

The width of the high-frequency gaussian is given by (2//i^)^/^. We require h'^ > 1 for proper 
application of steepest descent. This condition breaks down at lower frequencies, but we avoid this 
problem by combining the high-frequency results with the low and intermediate ones. We also note 
that one can show from eqn. (29) that our high frequency expression has a "super- Wien" behavior, 
i.e., it falls off faster than exp{—hi'/kT). 



5. Joint Spectrum 



To construct a joint spectrum that is valid in all three frequency regimes, we use the shell 
model. In order to recover the limiting cases for the shell radius and thickness, we sum the results 
from the low-intermediate and high frequencies. 



fm = Min < 



ukx 



l/(fcp+fcT-l) 



, Rc 



(32) 



rCexp(f>f^^i„^_i„t) 



(3 - fcp - kT)i'T(:Rc^^~^''~^^^ + kTU^^-k''^I^T 



~2 feo 
m, low— int 



<0k( V:i)^ c^i-^p 

^ ^'^'^ I, kp+l ) '^'m.high 



f2 
m, high 



Rc''" ') 



(33) 
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How do these results behave in the hmiting cases of very extended {Rc oo) and very compact 
(Rc — ^ 0) envelopes? For extended envelopes, the break frequency, inbreak ^ 0, so that the low- 
frequency regime disappears. Correspondingly, the shell radius at low and intermediate frequencies 
approaches the intermediate value, fm, iat (C/z/)^/^^. The high-frequency emission changes in 
this limit only insofar as Ti, increases somewhat. 

In the opposite limit of large optical depths, {Rch cxD) so that i?c — > 0), the photosphere 
approaches the cloud surface and the slope of the temperature profile becomes very steep {kr ^ 1). 
This case is difficult to realize in practice since we have assumed that there is no emission beyond 
Rc, and it is difficult, though not impossible, to realize such sharp boundaries in practice. A possible 
example of such a sharp boundary is a photoevap orating globule (Bertoldi &; McKcc 1990). We do 
not attempt to treat the Rc ^ limit in this paper, but instead confine our attention to i?c ^ 1- 

Figure 3a depicts the analytic SED (crosses) overplotted on the numerical SED (solid line) 
produced from DUSTY, with the three frequency regimes marked on the plot. This is a typical 
SED from the region of the physical parameter space labeled as "high mass protostar" in Figure 1, 
with L/M ~ 400Lq/Mq and S ~ 1 g cm~^. Figure 3b shows the combined contribution function. 
Figure 3c is a plot of the opacity curve, WDOl's Ry = 5.5. One should read these three plots left 
to right, i.e., follow the marked regions in the SED plot in Figure 3a and correlate them with the 
marked regions in the contribution function in Figure 3b. The spectral features in the contribution 
function in Figure 3b correlate with the spectral features in the opacity curve as depicted in Figure 
3c. For example, the 10 pLui (3 x 10^^ Hz) increase in the opacity translates to a corresponding 
increase in r^: as the r = 1 surface at this frequency is driven outwards, while the 5 /xm (6 x 10^^ 
Hz) decrease in the opacity causes to move inwards. 

6. Temperature Profile 

Numerical radiative transfer schemes solve for the temperature profile within the envelope by 
enforcing the condition of radiative equilibrium within every resolution element, i.e., that the total 
energy absorbed by a differential volume element equal the total energy emitted. This condition 
of energy balance is equivalent to the condition of zero flux divergence when the radiation field is 
time- independent (Mihalas, p. 48). In particular, in spherical geometry, this condition is equivalent 
to the constancy of luminosity as a function of radius. Approximate numerical radiative transfer 
schemes, such as that of Adams & Shu (1985), have enforced radiative equilibrium at a discrete 
number of points, using a particular form of the temperature profile, and shown that the total 
luminosity transported is approximately constant as a function of radius. Subsequently, more 
precise numerical radiative transfer solutions have solved for the temperature variation in the 
envelope in full generality, and enforced radiative equilibrium at fine resolution intervals, thereby 
guaranteeing virtually exact constancy of luminosity (IE97) . 

Here, in our analytic treatment of the radiative transfer problem, we have found that the far- 



IR emission can be inferred with good accuracy from a single power-law temperature profile. We 
determine the slope of this temperature profile from the self-consistency condition that the input 
luminosity exactly equal the emergent luminosity: 



(34) 



Inserting the shell expression for (eq. 14) into this expression gives 



QOjkp - 1) 
7r4 



/ 



Kiyi>^fm "Af^exp - 



kT{rm) 



- Ty{rm) di> . 



(35) 



This is analogous to requiring constancy of luminosity at two discrete spatial points; numerical 
schemes achieve greater accuracy by iteratively imposing zero flux divergence over fine resolution 
elements. Since and Arm depend on both kx and C, this equation gives a zeroth-order, self- 
consistent determination of the slope of the temperature profile, kx, when the fitting parameters 
L and C are taken to be equal to unity. To secure better agreement with the numerical results 
from DUSTY, we have simultaneously solved for kr, L, and C by imposing (35), and required that 
the maximum difference between the analytic and numerical SEDs between 3 mm and 30 /xm be 
as small as possible. Two of these parameters {kr and L) characterize the temperature profile and 
one (C) determines the break frequency separating low and intermediate frequencies. We find that 
the parameters L and C are indeed of order unity. 

For a given dust opacity, the shape of the SED, and therefore the values of the three fitting 
parameters, depend on three source parameters, Rc, T^i and kp. The dependence on T^i arises from 
the opacity, k.,^. We now recall that the opacity is scale- free longwards of 30 fim. If most of the flux 
is emitted longwards of 30 fim. (corresponding to hv/k = ASOK), then (x = from equation 
(15). In order to have most of the flux emitted longwards of 30 /xm, we require T^h < 300 K 
(for kp = 1.1, the condition is more stringent, Teh < 250 K). For lower values of T^h, the power 
law approximation to the opacity is valid and our fitting parameters will be independent of Teh- 
By comparing with DUSTY, we have found that it is possible to further simplify the functional 
dependences to 



With these dependences, we find an excellent level of agreement between our analytic results 
and the numerical results from DUSTY over a parameter space that spans low-mass protostars to 
ULIRGs. 

How do we expect kx and L to depend on i^c? We may heuristically understand the functional 
form of kx by considering the equation of radiative equilibrium, which requires that dust grains 
radiate as much as they absorb: 



kr = kriRc, kp), L = L{Rc), and C = C{kp) . 



(36) 




(37) 



where Bi, is the Planck function and Ji, is the mean intensity of the radiation field. 



-18- 



First consider the outer envelope, assumed to be optically thin. Insofar as we can make the 
approximation k^, oc in the LHS of this equation, it scales as T^+P-^ the RHS scales as 
in the limit of negligible optical depths. Therefore, in the outer parts of extended, optically thin 
envelopes, we have T{r) cx r~'^/^'^'^^\ If the envelope is cool (Tch ^ 100 K), we have p ~ 2 so that 
the slope of the temperature profile — > |. In fact, extended envelopes, which have large values of Rc, 
generally have higher temperatures, so that the effective value of p is somewhat less than 2 and the 
slope is somewhat greater than |. Next consider the inner, opaque region of the envelope. There 
the diffusion approximation is appropriate, and it gives T{r) oc r~^^^''>'^^^^^f^\ which is steeper 
than the slope in the outer envelope. For very opaque envelopes, the gradient can become steeper 
than this near the photosphere, just as in the case of a stellar atmosphere. We wish to choose a 
single power law, kr, to represent the temperature gradient near and outside the photosphere. For 
extended envelopes (large Rc) we expect kr to approach 2/(4 + whereas for compact envelopes 
(small Rc) we expect kr to increase as Rc decreases. We find that this expected variation can be 
represented by the sum of two inverse power laws, 

A R 

kT = ^ + ^ , (38) 

tic -n-c 

where ni and n2 are positive numbers, with ni being a small power in order to represent the 
dependence in extended envelopes. 

The normalization of the temperature profile is regulated by L, since 

(fcp-i) 

TchOcL 2/3+4(fcp-i) (39) 

from equation (7). This dependence is weak, but enables us to improve the accuracy of our fits - 
analogous to a temperature correction procedure. We may understand the variation of L with Rc 
as describing the transition from a modified blackbody to a protostellar envelope: In the limit of 
large optical depths, (small Rc) L is about 1 since most of the emission comes from the vicinity 
of a Rosseland photosphere. On the other hand, as Rc tends to infinity, the photosphere is less 
well-defined. In this limit, L is larger than unity, reflecting the effective increase in the total 
emitting area, as the intermediate frequencies probe the source function over the extended envelope. 
We depict the contribution functions for these two limiting cases in Figure 4. For large optical 
depths, fm is approximately constant, since most of the emission comes from the surface and the 
photosphere is relatively well defined. On the other hand, at large Rc there is a wide range of 
in the intermediate frequency regime. We find that this behavior can be represented by 

L = AlK , (40) 

with sufficient accuracy. 

To determine the values of kT, L, and C for given values of Rc and kp, we use the downhill 
simplex method of Nelder and Mead (1965) to minimize the largest error between the shape of 
the normalized analytic spectrum and the normalized DUSTY spectrum. We ran DUSTY at an 
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intermediate value of L/M (= AQLq/Mq), for kp = 1.1, 1.5, and 2, increasing E until we reached 
= 200K. To avoid calibrating our functional forms of k^, L, and C with respect to higher 
values of Teh, we then moved down the T^i = 200K isotherm towards large Rc until we reached 
S = 0.03. Along this line of constant L/M and this isotherm, we have ensured that the input 
luminosity equals the emergent luminosity. This is not exactly true away from this trajectory due 
to the weak dependence of the parameters on Teh that we have ignored. 

The temperature profile is determined by kx and L. Our results for the slope of the temperature 
profile are: 

0.48A:0-°5 O.lfc^-5 

i2e " Rc ' 

Figure 5 depicts fcr as a function of Rc- For L we find that 

L = 0.87i?°-°^^ (42) 

is sufficiently accurate for all kp. 

We find that setting C = constant for each value of kp allows for sufficient accuracy (in fact, 
letting C vary as a function of Rc does not increase accuracy). Values of C{kp) are selected to give 
good agreement with the break frequency for large R^ envelopes. This gives C = 1,0.9, and 0.5 for 
kp = 1.1, 1.5 and 2 respectively. For other values of kp in the range 1 < fcp < 2, C can be found 
from 

C = 0.27 + l.Zkp - OMp"^. (43) 

A more detailed discussion of this factor and its determination from the spectrum itself, i.e., the 
break frequency, for extended envelopes is given in §8. 

We conclude with a caveat on our temperature specification procedure. The near-IR flux 
depends sensitively on the temperature profile, and we cannot recover the near-IR flux accurately 
with a single power law for the temperature profile. For example, as shown in Figure 3, it is clear 
that the emission at 5//m, where there is an opacity minimum, originates inside the characteristic 
radius. Our single power law underestimates the temperature, and therefore the emission, for 
regions well inwards of the photosphere. This problem could be remedied with a two-component 
power law for the temperature profile, but that is not necessary here since we are focusing on 
emission at longer than near-IR wavelengths. 



7. The Analytic SED and Its Accuracy 



For convenience, we collect the equations that we have derived for the far-IR SED of dusty 
sources: The mm to far-IR SED is given by: 



= 16Tr'^{kp - 1)R, 



ch 



2hu^ 



hv 



kT{fm) 



(44) 
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The characteristic emission radius, fm = fm{^), is the location in the shell where most of the flux 
in a given frequency-band originates from, (the "m" is for maximum), and is given by: 



Min (f 



m,low— int '"m,high 



(45) 



where the total Vm is the sum of the high frequency and the combined low-intermediate frequency 



''m,high 



''mjlow— int 



Ri}l/kT + cy^T 



(46) 
(47) 



The shell thickness, A.rm is given by: 



Arr. 



{S-kp- kT)i>TCR, 
■i 



-(S-kp-kr) 



f'^ kp 

m, low— int 



1 -I- AKi, y j I'm, high 'm,high-"'C ; 



(48) 



where the argument {2> — kp)/kT for the Gamma and Zeta functions has been suppressed for clarity, 
and h'l^ is given by Eqn. 30. The ratio of hv to kT equals C for the intermediate frequencies and 
is given by: 



C = 0.26 + 1.3A;. - 0.6A;^ 



and 



kx = 



0.48fc0' 



0.1fc5-5 

~0.02fci09 + p0.7fci-9 
Kc ric 



(49) 
(50) 



The SED variables and the source parameters are related by: 



Rr. 



{ 



4+/3 



{L/M)^~ 
ia X 0.87 



(3 kp)Ku^^ 



4(fcp - 1)^0^ 



■ 1.92/3+4(fc|0-l) 



(51) 



^ch 



L/M 



2.92-fep 

4(7 X 0.87S"^^ 



4(fcp - 1)^0^ 
(3 — kp)Hi,Q 



4{fcp-l)+1.92/3 



where we have substituted for the values Al and a from (42). 
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The physical parameter space in which we have determined the accuracy of our analytic SED, 
as compared with the numerical solutions from DUSTY, is given in Figures 1, 6, and 7 for three 
density profiles, kp = 1.1, 1.5 and 2. The parameter space extends over luminosity-to-mass ratios 
L/M = {0.1 — 4000) Lq/Mq and surface densities E = 0.01 — 100 g cm~^, which encompass the full 



- 21 - 



range of astronomical dusty sources. We restrict our attention to T^h < 300 K (250 K for kp = 1.1) 

so that most of the emission is longwards of 30 ^m and the slope of the temperature profile is 
approximately independent of Tch- We also restrict our attention to i?c > 1, since our assumption 
of a spherical cloud with a sharp boundary is likely to break down for smaller envelopes. Within 
these boundaries, the SED given by the above equations is accurate to within a factor 2 at worst. 
The largest errors occur for kp = 1.1; for kp = 1.5, the SED is accurate to within a factor of 1.6, 
and for kp = 2 it is accurate to within a factor 1.5. If we focus on surface densities in the range 
0.01 — 3 g cm~^, which is indeed where the majority of astrophysical source lie, the accuracy is 1.5 
for kp = 1.5 and 1.3 for kp = 2. We also note that our equations hold for kp = 1, as numerical 
runs verify that kp = 1 envelopes are well approximated by /cp = 1.1 (to better than 10 %) over our 
parameter space. Finally, we note that our largest errors, as compared to the numerical results, 
are at 30 /xm, on the low Rc end (large S region), where there is very little flux. As a result, our 
accuracy in inferring source parameters (given that our assumptions stated in §2 arc satisfied) is 
generally better than a factor of 1.5, a point that discuss in more detail in Paper II. SED fits, to 
a low-mass protostar, massive protostar, and ULIRG, respectively, are given in Figure 2a, Figure 
2b, and Figure 2c. 

We note here briefly (see Paper II for a more detailed explanation of the application of our 
methodology to observed data) that our accuracy criteria hold only so long as our assumptions are 
valid. In particular, many of the sources we consider in Paper II have significant high frequency 
fluxes that our solution cannot account for. High frequency emission can escape due to inhomo- 
geneities in the envelope, due to the presence of an accretion disk, or from distributed sources of 
luminosity in the field of view that are not obscured by dust. We consider the effects of an inhomo- 
gcncous dust distribution in Paper III, and the effects of distributed sources of luminosity in Paper 
IV. To avoid fitting to data that may predominantly arise from one of these additional agents, in 
Paper II we have performed our fits from mm - 60 fxm, where these additional complexities are 
unlikely to significantly change the emitted spectrum. If the fit based on the mm - 60 /xm data also 
fits the high frequency data, that is a good indication that the high frequency data are not due to 
these additional agents. 

Our basic assumption of spherical symmetry can be violated by the presence of a disk that is 
large and massive enough to influence the far-IR spectrum. Since the disk radius approximately de- 
marcates the region where the density, and therefore the temperature, profiles become nonspherical, 
our method is applicable only to the case in which the disk radius is smaller than the characteristic 
radius, < Rch- Even if this geometric condition is satisfied, we also require that emission at 
long wavelengths (where the emissivity scales as the product of the mass times the temperature) 
from the disk be small compared to that from the envelope, which implies 



We have also assumed that emission from the vicinity of the dust destruction radius at R^d 
does not significantly influence the far-IR spectrum. A sufl&cient condition to ensure this is that the 
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optical depth to the dust destruction front at photospheric frequencies be large; we require r > 3 
at the frequency t'peak at which uLn peaks, which implies 

< ^ T- ■ (54) 



1 + 3 Vak^c" 



Figure 8 depicts Upeak = ^peak/i^ch as a function of Rc, and §8 gives corresponding fits to the 
curves. This optical depth condition is satisfied for kp = 2 and 3/2 for all Rc in our parameter 
space. It is also satisfied for kp = 1.1 when < 10. However, for kp = 1.1 and Rc ^ 10, the 
optical depth to the dust destruction front at t'pcak is not large. Nevertheless, emission from i?dd 
does not influence the spectrum because the far-IR luminosity originating from the vicinity of the 
dust destruction radius is small compared to that from the vicinity of the characteristic radius in 
our parameter space. An upper limit on emission from the vicinity of the dust destruction front is 
a blackbody radiating at Tdd- Hence, at intermediate frequencies, a sufficient condition for ignoring 
emission from the dust destruction front is 

TadRld « TchRch exp (-C) , (55) 

(with hu/kT = C), which is always satisfied in our parameter space. At frequencies well above 
i^peak, where the emission from the vicinity of the dust destruction front is significant, the optical 
depth is large enough to suppress it. 



8. Shape of SED: Rc and T^h 

As long as the assumptions stated in §2 & §7 are satisfied, i.e., given a spherically symmetric, 
homogeneous distribution of dust illuminated by a central source of radiation, such that the emitted 
spectrum is not significantly affected by emission from the dust destruction front, the shape of the 
long-wavclcngth SED will depend only on Rf. for a given density profile and dust model. Thus, 
moving along a line of constant Rc preserves the shape of the SED, while shifting the peak of the 
SED as one intersects lines of different Teh. The ratio of the peak frequency to the characteristic 
frequency is shown in Figure 8 and can be approximated as 

!^~0 82fc I 5-4-1.8fcp 

^ch. He 

Similarly, moving along an isotherm, i.e., a line of constant Tch, changes the shape of the SED as 
Rc changes but preserves near constancy of the peak of the SED since the variation of fpeak with 
Rc is weak over most of the parameter space (see eq. 56). 

Note that the ratio of the peak frequency to the characteristic frequency decreases as Rc 
increases because the intermediate frequencies, which are emitted at lower temperatures, become 
relatively more important in the total energy balance. In the low Rc regime (large optical depths). 
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our results are qualitatively similar to the blackbody limit, as z^peak ~ ^^ch (recall that Vpeak is the 
frequency at which lyL^^ peaks). Over much of the astrophysical parameter space, however, spectra 
peak at fpeak ~ (1 - 1.6)fch- 

Low mass protostars, as shown in Figure 9a, are characterized by large i?,,. and have broad 
SEDs. They have extended envelopes with temperatures low enough such that the intermediate 
frequency region, with its flatter slope, shows a clean separation from the low-frequency region. On 
the other hand, the far-IR spectra of massive protostars, ULIRGs and super-star clusters resemble 
quasi-blackbody spectra, where the low, intermediate, and high-frequency components have become 
smeared together (Figure 9b). Fits to observed data (Figure 2a, Figure 2b, and 2c) also display 
this variation of shape, as Rc varies over the parameter space, from the low-mass protostars to the 
massive protostars. 

We also give the general relation between the characteristic parameters and the photospheric 
parameters. As noted earlier, our definition of the characteristic parameters coincides with photo- 
spheric parameters when core radius is much larger than i?ch- However, as Rc becomes small, the 
Rosseland photosphere moves to the surface of the core. Therefore, we approximate the relation 
between the photospheric and chracteristic parameters as a harmonic mean of Rc and Ret- 

R..-^. (57) 

We use L = 47ri?pjja"rpj^ = AnR^^aT^^^L to see that the photospheric temperature varies as: 

rph = rehLV4(i + ^-1)1/2. (58) 



8.1. The Three-Component Spectrum of Extended Envelopes 

We now consider the limiting case of highly extended envelopes, i.e., Rc oo. We reach this 
region of astrophysical parameter space by moving downward and to the left along an isotherm 
towards low L/M and low S, i.e., past low mass protostars towards very low luminosity cores, or 
nearly prestellar cores. Figure 11 depicts a large Rc envelope and its corresponding contribution 
function, with the break frequency marked in both diagrams. We note that in this limiting case, 
since kr — > constant, we may infer the source parameters from fpeak s^nd i^reak only. The peak 
frequency is given by equation 56; the break frequency is given in terms of L/M and S by 



inbreak — C < 



(L/M)S(4+/3)//3 



4C70.87 



(3 



,4(fe,-TO 



4//3- 



/3fcT/[1.92^+4(fep-l)] 



(59) 



where C is given in equation 43. Figure 10 depicts the factor C{Rc) for the three standard density 
profiles; as shown, C — > constant for large Rc- 
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The clean separation of the intermediate frequency regime in this limiting case allows us to 
specify the temperature dependence of the spectrum in the intermediate frequency regime, 

oc T^-^''" , (60) 

which follows from Eqn. 12 and Eqn. 17 with fc^- = 1/3 and (3 = 2. The slope of kp = 1.1, 1.5 and 2 
envelopes depends on T^^, T^/^, and Teh respectively, so that fcp = 1.1 envelopes are most sensitive 
to the temperature. 

As we discuss in §8.3, while it is not generally feasible to discriminate density profiles from 
the far-IR spectra of unresolved sources, it is possible to do so for extended envelopes. In this 
case, we may write down a simple expression for the intermediate frequency slope, which holds for 

oc j.3+/3-(3-fe,)A^(Re) _ (61) 

For envelopes that are sufficiently extended {Rc > 2000), one can infer kp directly from the measured 
slope in the intermediate frequency regime. The pristine environments of isolated, low-luminosity 
cores allow for the most robust determination of the density profile from the far-IR SED, without 
requiring resolved observations of the source itself. Such a determination of the density profiles 
for these sources would also provide a direct observational test of the star formation scenario, as 
different theories of star formation are characterized by different density profiles. 



8.2. Density Profile 

The inference of density profiles from the far-IR SED has been treated in various ways in 
the literature. Some authors have determined the density profile from fitting to the SED and 
the intensity profile concomitantly (Shirley et al 2002, Mueller ct al 2002, Beuther et al 2002), 
which works if the intensity is solved for self-consistently. The Adams (1991) approximation for 
the spatial distribution of the millimeter and submillimeter emission, and its use to infer density 
profiles (e.g. Beuther et al 2002), is based on the assumption that l/3</cr^2/5, which is a good 
approximation for large Rc- However, the dense envelopes of massive protostars imply they are 
characterized by lower values of Rc, and therefore higher values of kx- Hence, the accuracy of the 
Adams approximation for dense envelopes is degraded in this region of the parameter space. Some 
authors (e.g. van dcr Tak ct al 2000, OLD99) have claimed that density profiles can be inferred 
from the observed SED alone. Others, (e.g. Correia et al 2004) have noted that varying density 
profiles give equally good fits to the SED. Therefore, the question remains: under what conditions 
can the density profile be inferred from the mm to far-IR SED alone? 

So far, we have shown that for a given density profile and dust model, the SED is a function 
only of L/M and S. To address the question of whether inference of density profiles is feasible, we 
now ask under what conditions can one distinguish the predicted mm to far-IR fluxes for envelopes 
characterized by different density profiles? To do this, we consider SEDs produced by envelopes 
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characterized by the same L/M and E, but with different density profiles. Figure 12a depicts the 

the 1mm to 60;um flux ratio vs. the the 30;um to GO^um flux ratio for the three density profiles, 
kp = 1.1,3/2,2, with each trio having the same L/M and S values. The kp = 3/2 points (crosses) 
correspond to an isotherm, Teh = 210 K, with the bottommost point corresponding to Rc = 43 
and the topmost point to Rc = 370. Note that for the first two trios, the colors for the different 
density profiles are not clearly distinguishable (i.e., they do not differ by more than a a factor of 
two), while for the third trio (Rc = 180, Teh = 210 K for kp = 3/2) they arc. Thus, the different 
density profiles are distinguishable at sufficiently large Rc. The region in the L/M — S plane in 
which at least one of the two colors differs by at least a factor 2 between the kp = 1.1 and kp = 2 
cases is shown in Figure 12b. 

The basic intuition behind this result is that envelopes with large values of Rc (i.e., envelopes 
in which the photosphere is small compared to the size of the core) have an extended range of 
intermediate frequency emission in the far-IR. Emission at intermediate frequencies originates from 
a range of radii and temperatures, and as a result, the emission from large Rc envelopes depends 
more sensitively on the density profile than that from low or moderate Rc envelopes. This is true 
for any value of T^i, i.e., it will be easier to discriminate density profiles for larger values of Rc 
whether we consider sources at high T^ii or low T^h- However, as shown in Figure 12b, there is 
some dependence on T^h, since the ratio of the 30 /xm emission to 60 //m emission is sensitive to the 
temperature; this effect is ameliorated by the fact i^peak/^^ch increases at low values of Rc, where 
Tf-h is low also. 

We conclude that while it may be feasible to infer the density profile from the SED for low 
and intermediate mass protostellar sources, it is difficult to do so for high-mass protostars and 
extragalactic sources as they are generally on the low Rc end of the parameter space. 

9. Conclusion 

We have presented an analytic, self-consistent solution for the spectral energy distribution 
of homogenous, spherically symmetric, dust-enshrouded central sources of radiation that are not 
affected by emission from the dust destruction radius. The main points are: 

1. For a given dust model and density profile, the SED is determined by three parameters, 
the luminosity of the central source, L, the mass of the envelope, AI, and the size of the enve- 
lope, Rc- The shape of the SED is distance independent, and is determined by the two distance- 
independent parameters L/M and S = M/{nR^). The cases of greatest relevance in astrophysics 
have O.IL0/M0 < L/M < 4OOOL0/M0 and 0.01 g cm'^ < S < 100 g cm^^. We consider power- 
law density profiles, p oc r~^p, with 1 <kp <2. 

2. The characteristic radius, i?ch) is analogous to a Rosseland photosphere. The characteristic 
temperature, Tph, is given by L = A^irLR^^aT^^, where L is a number of order unity. The emission 
at a frequency v can be viewed as originating from a shell centered at rm{v), with thickness Ar^(z/) 
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and with an attenuation e Most of the emission originates outside Rch, since for hv > kT^.^ the 
optical depth is large, whereas for hu < kT^h most of the emitting mass is outside the photosphere. 

3. The SED is readily described in terms of the characteristic temperature, Tph, and the 
dimensionlcss parameter Rc = Rc/Rch- For a given dust model and density profile, T^h and Rc are 
unique functions of L/M and S. Low-mass protostars have extended envelopes with core radii that 
are large compared to the photospheric radius, so that i?c S> 1. High-mass protostars, ULIRGs 
and super-star clusters generally have less extended envelopes (smaller values of Rc)- The width of 
the SED correlates with Rc'. sources with small values of Rc have quasi-blackbody SEDs, whereas 
those with large values of Rc have broad SEDs. 

4. The SED has three frequency regimes. Low frequencies are emitted by the entire envelope, 
and have hv < kT{Rc). For the density profiles we have considered, most of the low-frequency 
emission comes from the outer part of the envelope; it is optically thin. Intermediate frequencies 
have kT{Rc) <hu < kT^i, so that they are suppressed in the outer envelope; most of the emission 
is optically thin. Finally, high frequencies have hv > kT^h- The flux at high frequencies is the result 
of a competition between the opacity, which favors emission from larger radii, and the temperature, 
which favors emission from smaller radii, where the dust is warmer. Low frequencies vary as 
intermediate frequencies vary as i'^+f^~i^-'^p)/''T{Rc) ^ and high frequencies have a "super-Wien" 
behavior - they fall off faster than e^p{—hu/kT). 

5. We approximated the temperature profile near the photosphere as a power law, T oc r~'^'^, 
where the temperature gradient, fc^, is determined from imposing the self-consistency condition 
that the input luminosity exactly equal the emergent luminosity. Wc then showed that adoption 
of this approximate temperature profile leads to good agreement (usually within a factor ~ 1.5 
between 3 mm and 30 iim) with spectra calculated with the numerical radiative transfer code 
DUSTY. 

6. Longwards of 30 /nm, where the opacity is scale free, the shape of the SED is a function only 
of Rc and kp. However, the dependence on the density profile (kp) is significant only for relatively 
large values of Rc- Large Rc envelopes, such as those of low-luminosity, low-mass protostars, present 
an opportunity to infer the source parameters from two observable frequencies, fpeak and i^reak) 
where t'brcak divides the low and intermediate frequency regimes. At large Rc, the intermediate 
frequency regime has a slope that depends explicitly on kp. 

We thank Bruce Draine, Eugene Chiang, David HoUenbach, Xander Tielens, and especially 
Gibor Basri for helpful discussions. We also thank Maia Nenkova and Dejan Vinkovic for answering 
innumerable DUSTY-related questions. This research is supported in part by NSF grant ASTOO- 
98365. CFM gratefully acknowledges the support of the Miller Foundation for Basic Research. SC 
gratefully acknowledges a dissertation year fellowship from the Graduate Opportunity Program. 
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A. Relation to DUSTY Parameters 

We give the relations between the characteristic parameters and the DUSTY parameters. 
First, we note that the DUSTY dimensionless parameters are: y = Rc/Raa = Rc/ Rdd', Ty^j the 
optical depth through the shell at frequency vq; and kp, the slope of the density profile. The 
dimensional parameters are the dust dcstrTiction temperature, Tdd, and the temperature of the 
central illuminating source, T*. In addition, the shape of the opacity curve is required as input. 
The dust destruction radius, R^Ai can be expressed as 



R, 



■dd 



(Al) 



as in IE97. 



Thus, the value of the dust destruction radius depends explicitly on the magnitude of the 
luminosity. For our purposes, the parameters Tdd and i?dd do not enter into the problem, as we 
consider optical depths to the dust destruction radius large enough such that it is opaque at far-IR 
wavelengths. We choose vq to be in the power-law portion of the opacity curve so that we may 
express Tj/q as 



where Eaa = i?dd/^ch- 



B. Escape Fraction 

Here, we give the derivation of the escape fraction for spherical geometry: 

/(r)= r , (Bl) 

where ^ = cos{9). For escape at small angles (which is appropriate for large optical depth), we 
evaluate t{0) for ^ <C 1. We express the optical depth as: 

Ty = h^ikp -I) j r'-^^dS , (B2) 

where the path length is equal to 



with r as the source point. We evaluate the optical depth for ^ ^ 1, to find that: 
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Now, we evaluate the integral for the escape fraction with respect to angle and find that: 

f{r) = ^ ^ r,(r,0)»l. (B5) 

We assume that the Tjy(r, 0) = 1 surface is at r ^ Re, so that the escape fraction is unity near the 
surface of the core. We recover the proper limit when r = Rc, i.e., we recover / = 1 in that case, 
if we modify this to be: 

e-T-,.(r,o) 
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Fig. 1. — {a)L/M[L(7)/MQ] vs S[g/cin^] for fiducial density profile, kp = 3/2. Teh is the effective 
photospheric temperature and Rc is the ratio of the core radius to Rch, which is like the Rosseland 
photosphere (b) with change in opacity normalization of a factor of two (dashed line), due to 
ice-mantles 
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IRAS 16293-2422 , IRAS 18454-0136 ^ UGC 5101 




Fig. 2.— (a) SED fit to low-mass protostar, L/M ~ 7Lq/Mq, S ~ 0.3g cm-^ (b) SED fit to 
massive protostar, L/M ~ IILq/Mq, S ~ 1.5g cm^^ (c) SED fit to ULIRG, L/M ~ SSLq/M©, 
S ~ 0.3g cm~^. Details of fitting procedure and inferred parameters given in Paper II 
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Fig. 3. — (a) SED for typical high mass protostar with frequency regimes marked (sohd Une is 
DUSTY SED and crosses analytic SED) (b) Contribution function (the characteristic emission 
radius), in dimensionlcss units, f^, with frequency regimes marked, (c) WDOl opacity curve. The 
spectral features in the SED and opacity curve as shown in (a) and (c), e.g. the 3 x 10^'^Hz (10//m) 
absorption feature, correlate with the location in the envelope this emission is coming from, as 
shown in (b) 



-34- 



1000.0 =- 



100.0 - 



10.0 



1 .0 ^ 




0.1 

10 



1 1 



10 



12 



10 



13 



'[Hz] 



Fig. 4. — Contribution function (the characteristic emission radius), in dimensionless units, f^, for 
a compact envelope, Rc = 4: (solid line), and for an extended envelope, Rc = 2000 (dotted line) 
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Fig. 5. — (a)Slope of temperature profile, fcy as a function of Rc- Asterisks denote kp = 1, crosses 
denote kp = 3/2, and diamonds denote kp = 2 
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Fig. 6.— L/M[Lq/Mq] vs E[g/cm2] for kp = 2 
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Fig. 7.— L/M[Lq/Mq\ vs E[g/cm2] for kp = 1 
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Fig. 8. — The ratio of the peak frequency to the characteristic frequency, fpeak/^ch^ ^ a function 
of Rc- Asterisks denote kp = l, crosses denote A;^ = 3/2, and diamonds denote kp = 2. 
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Fig. 10. — C{Rc) for kp = 3/2 (crosses), kp = 1 (asterisks), and kp = 2 (diamonds) 
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Fig. 11. — (a) SED for large Rc envelope (extended envelope) with break frequency marked, (b) 
Contribution function 
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Fig. 12. — (a)F(lmm)/F(60/im) vs F(30/im)/F(60/Ltm) for kp = 2 (diamonds), kp = 3/2 (crosses), 
kp = 1 (asterisks); kp = 3/2 curve corresponds to Teh = 210 isotherm, (b) L/M, S plot for kp = 3/2, 
dashed Une demarcates region below which density profiles can be inferred from far-IR SED 



